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Effects Of Granulation Upon Larger-Scale Convection 
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Abstract. We examine the role of small-scale granulation in helping to drive super- 
granulation and even larger scales of convection. The granulation is modeled as localized 
cooling events introduced at the upper boundary of a 3-D simulation of compressible 
convection in a rotating spherical shell segment. With a sufficient number of stochas- 
tic cooling events compared to uniform cooling, we find that supergranular scales are 
realized, along with a differential rotation that becomes increasingly solar-like. 



1. Introduction 

The solar differential rotation profile e xhibits prominent rad ial shear layers near the 
top and bottom of the convection zone (|Thompson et alJll99d ). The near-surface shear 



layer (NSSL) occupies the upper 5% of the Sun by radius, whereas the tachocline be- 
gins near the base of the convection zone. Much of the dynamics of the NSSL is due 
to the vigorous granular-scale convection that is driven by radiative cooling and large 
superadiabatic gradients. The collective interaction of these granular-scale flows (av- 
erage sizes of 1 Mm, lifetimes of 0.2 hr) is a major component in the formati on of 
supe r granular (15-35 Mm , 24 hr) and mesogranular (5-10 Mm, 5 hr) scales (e.g. iRast 



20031 : Asplund et al. I [20091 ). The preferential influence of the Sun's rotation upon su- 



pergranular and larger scales may lead to the anisotropic angular momentum transfer 
seen in the NSSL. The overturning times and sizes of supergranules suggest that gran- 
ular convection is weakly influenced by rotational effects. The interactions between the 
small-scale convection, larger velocity gradients, and rotation yield Re ynolds stresses 



which are significant in transporting angular momentum within the layer (jDeRosa et al 
20021 1 



A direct numerical simulation of all these dynamical scales is not possible given 
current computational resources. Therefore, in order to investigate the dynamics of 
the 100-1000s of supergranules necessary for adequate statistics, we must model the 
small scales that we cannot capture. We examine two methods of carrying the required 
energy flux through the top of the domain, one using a spherically symmetric (uniform) 
cooling function while the other is a granular-scale stochastic cooling model. In this 
paper, we investigate the dynamics and differential rotation profiles achieved in our 
simulations. The governing equations and numerical approach used in solving them are 
briefly discussed in £}2j Our granulation model is reviewed in £}3j while the connection 
between these simulations and the near-surface shear layer of the solar convection zone 
is discussed in 



2. Formulating the Problem 

As a compromise between physical accuracy and more rapid time evolution, a perfect 
gas and a radiative diffusion approximation are employed in our simulations. The spher- 
ical segment domains used in our simulations involve large portions of the Sun's inherent 
spherical geometry, which is necessary to properly capture the effects of rotation on 
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supergranular scales. To model these dynamics below the solar photosphere, we use a 
heavily modified version of a fully compressible 3-D hydrodynamic simulation code. The 
Curved Shell Segment (CSS) code is a mature modeling tool whi ch solves the Navier- 
Stoke s equations of motion in rotating spherical segments (e.g. iHurlburt &: DeRosa 
2008). To simulate the larger scales of motion that are likely to occur in the solar 
convection zone, a large-eddy simulation (LES) model is employed. The scales that 
are not explicitly computed in these simulations are parametrized and included in a 
sub-grid scale model of turbulent transport and diffusion. The equations of motion 
solved in CSS are: 
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where the symbols p, u, P, T, S, and V are the density, velocity, pressure, temperature, 
specific entropy, and viscous stress tensor respectively; the equation of state is P = 
ple S ' Cy ; f2 is the angular velocity of the rotating frame; ks is a turbulent-eddy entropy 
diffusion and K r is the thermal diffusion; g is the local acceleration due to gravity; eo is a 
source term that carries the required flux through the boundaries. Here p, k$, K r , g, and 



eo are functions of radius only. The diffusive energy flux is q 
The granular-scale cooling is contained in q^. 

2.1. Numerical Methods 
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The equations of motion (Eqn. [T]) are evolved on a uniform spatial mesh. Temporal dis- 
cretization is accomplished using an explicit fourth-order accurate Bulirsch-Stoer time- 
stepping scheme. The spatial de rivatives a re computed using a modified sixth-order 
compact finite difference scheme ( Leldfl992l ). A 3-D domain decomposition divides the 
full spatial mesh into sub-domains. The boundary information necessary to compute 
spatial derivatives is passed between nearest-neighbor sub-domains using MPI, while 
computations are shared among the master and slave cores within a supercomputer 
node using OpenMP. 

To begin to understand the dynamics of the near-surface shear layer, we have 
constructed numerical simulations in spherical shell segments that encompass most of 
the near-surface shear layer and some of the deep interior. A stellar evolution code 
is used to establish a realistic initial stratification for the simulations. Since a perfect 
gas is assumed, the He and H ionization zones can not currently be simulated. Thus, 
the upper boundary (rg) is taken to be 0.995 Rq in order to exclude these zones. 
The lower radial boundary (r\) is the radius which corresponds to a given density 
contrast (Ap = pij pi, where p\ = p{r\) and P2 = p{r2))- The momentum and entropy 
diffusivities (p and k s ) are calculated based upon the desired Rayleigh number (Ra) at 
r2 and so that the Prandtl number (Pr) is unity throughout the domain. We choose 
to model these sub-grid-scale diffusivities based upon mixing-length theory, so that the 
diffusivities become functions of radius. 



3. Modeling Granular Convection 

The strong downflows that occur at edges of a granule drag more dense, lower entropy 
fluid into the interior. Hence, in our simple model we assume that the net mass flux is 
zero and the net entropy flux is negative when integrated over the volume of a granule. 
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It is also assumed that the collective action of all the granular-scale dynamics contained 
within our domain transports the required energy flux through the upper boundary. 
Therefore, we have modeled the granules as cooling events that are advected with the 
flow, last for an average of 10 minutes, and have as small a spatial scale as the grid 
and diffusion will allow. The spatial structure of a cooling event is evident in Eqn. 
[5J The amplitudes of the cooling events are chosen at random from an exponential 
distribution with a mean determined by Eqn. [3j The horizontal extent of the cooling 
events is normally distributed with a mean angular size of 2A0 ((<r)) in 9 and <p. The 
spread of the radial exponential is fixed to ay = 2Ar to speed computation. The cooling 
function is thus 

N(t) 
£ Qr, 



-V • (q R r) 



+ 



1 



-{ri- 



(2) 



where Q n (t) is the amplitude of the cooling event at time t, and the cooling events are 



also horizontally advected with the flow: 6 n (t) 
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A mean number of cooling events ({N)), with ((cr)), determines the mean strength 
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Cooling events have a mean lifetime of 10 minutes, but can deviate significantly from 
this as the lifetimes are normally distributed as well. When a cooling event has reached 
the end of its lifetime, it is deleted and at least one new cooling event is generated. A 
uniform distribution is used to choose how many new events are generated. Similarly, 
at each time step there is also a chance to add new events. In future models, improved 
st atistics for th e gran ular scale cooling based upon the surface convection simulations 
of Tramp edach ( 20071 ) will be implemented. 



Table 1: Properties of the Computational Domain 





Case 1 


Case 2 


Case 3 


Case 4 


(N) 





250 


500 


1000 


(Ra) 


5.98 x 10 6 


7.63 x 10 6 


7.90 x 10 6 


7.94 x 10 6 


(Re) 


133 


152 


151 


150 


(Ro) 


13.3 


17.7 


16.9 


16.6 



Table 1. For each case, the depth of the convection zone as modeled is 60 Mm, Ap 
is 400, the angular size is 20° x 20°, with a resolution of 128 x 256 x 256. (N) is 
the mean number of cooling sites. All cases were run with the same diffusivities and 
initial background. (Ra), (Re), and (Ro), are the average Rayleigh, Reynolds, and 
Rossby numbers in the domain. 



4. Nature of the Convection 



In a typical simulation, the convective cells are broad upflows bordered by strong, 
narrow downflows as seen in Figure [Th. The upflows are warm and the downflows 
are cool relative to the mean temperature. The typical radial velocities in Figure [Th. 
are between —1400 ms _1 for downflows and 400 ms _1 for upflows. The spatial scales 
present in the four cases can be seen in power spectra (FigurefTJa). The supergranule-like 
scales in Case 1 are about 40 Mm, just outside the accepted range for supergranules. 
A spherically symmetric cooling is used in Case 1 which leads to horizontal scales 
near the surface that are primarily selected by the geometry of the domain. In cases 
with granular-scale cooling, the size scales seen in the simulation decrease due to the 
granular-scale forcing as evidenced by the leftward shift in the power of the peaks below 
40 Mm in Figure [TJa. In Cases 2-4, supergranular scales become much more prevalent 
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Figure 1. (a) A 20° x 20° snapshot of radial velocities (u r ) taken at 0.99 Rq in 
case 2: downflows are black, upflows are white, (b) Horizontal scale i in u r power 
spectra averaged in time and latitude from cases 4 (solid), 3 (dotted), 2 (dashed), 
and 1 (dash-dotted), (c) Horizontally averaged rotation rate fi(r), lines as above. 



with power increasing by 50-100% in the band between 10 Mm and 35 Mm. In Case 
2, a prominent peak occurs around 10 Mm due to the relative paucity of cooling sites. 
Since there are only 250 cooling sites on average, each site necessarily induces a strong 
downflow, which ensures that the enthalpy flux is sufficient to carry out the solar flux. 
Furthermore, there are typically 15 sites at a single longitude or latitude so that the 
average distance between them is 10-15 Mm. Thus, since each cooling site is likely to 
be a strong downflow, power is present on the 10 Mm scale. Such a peak does not occur 
to the same extent in Case 3 or Case 4 because the probability of inducing a strong 
downflow at a cooling site decreases by factor between 2 to 4. However, a larger (N) 
reduces the average distance between cooling sites and thus allows for more interaction 
of the flows generated at the granular scales due to the advection and coalescence of 
the cooling sites. 

As (N) is increased, the granular scale cooling disrupts the mean horizontal flows in 
the upper 20% of the domain which reduces the mean rotation rate as seen in Figure [1};. 
The reason the mean flows are affected is that when a cooling site forms in an upflow, 
it begins cooling nearby gas causing contracting horizontal flows and a downflow. So, 
the net effect of this cooling site is to reduce the horizontal flows within the upflow 
as well as slow the upflow's ascent as it has been cooled. Hence, as (N) increases 
the average number of cooling sites per upflow increases thus decreasing the mean 
horizontal flows as in Figure [1}:. The radial boundaries are currently impenetrable and 
stress-free leading to overturning of upflows at the top and downflows at the bottom of 
the domain. This is the primary reason the rotation rate at the bottom does not flatten 
out and why there is an unrealistic shear layer near the surface. Future simulations 
will have inflow-outflow radial boundaries which should eliminate these effects. 
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